<!DOCTYPE html>
<meta charset="utf-8">
<script type="text/javascript" src="http://cdn.mathjax.org/mathjax/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML">
</script>
<script>
  MathJax.Hub.Config({
                      tex2jax: {inlineMath: [['$', '$'], ['\\(','\\)']]},
                      TeX: { equationNumbers: {autoNumber: "AMS"} },
                      "HTML-CSS": { showMathMenu: false,
                                    scale: 90 }
                     });
</script>
<link rel="stylesheet" href="http://examples.gurobi.com/base.css">
<style>
a:hover.screenshot {
  opacity: .7;
}

.subunit-label {
  fill: #777;
  fill-opacity: .5;
  font-size: 20px;
  font-weight: 300;
  text-anchor: middle;
}

.place,
.place-label {
  fill: #444;
}

text {
  font-family: "Helvetica Neue", Helvetica, Arial, sans-serif;
  font-size: 12px;
  pointer-events: none;
}

</style>
<body>
  <ul id="nav">
    <li class="current"><a href="#intro">Intro</a></li>
    <li><a href="#problem">Problem</a></li>
    <li><a href="#model">Model</a></li>
    <li><a href="#implementation">Implementation</a></li>
    <li><a href="#demo">Live Demo</a></li>
    <li><a href="#try">Try Gurobi for Free</a></li>
  </ul>
  <div id="example_container">
    <div class="example_section" id="intro">
      <h1>The Traveling Salesman Problem</h1>
        <subtitle>with integer programming and Gurobi</subtitle>
    </div>

    <p>
      In this example we'll solve the Traveling Salesman Problem.
    </p>

    <p>
      We'll construct a mathematical model of the problem,
      implement this model in Gurobi's Python interface, and compute and
      visualize an optimal solution.
    </p>

    <p>
      Although your own business may not involve traveling salesmen, the
      same basic techniques used in this example can be used for many other
      applications like vehicle rounting, circuit design and DNA sequencing.
    </p>

    <h3>
      Click the screenshot to skip directly to the Live Demo!
    </h3>
    <p>
      <a href="#demo" class="screenshot">
        <img src="screenshot.png" alt="Live Demo" style="width: 100%; vertical-align: middle;">
      </a>
    </p>

    <div class="example_section" id="problem">
      <h2><a name="problem">Problem Description</a></h2>

      <p>
        The Traveling Salesman Problem (TSP) is a classic problem in
        combinatorial optimization. It was first formulated as an integer
        program by
        <a href = "http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.134.9319&rep=rep1&type=pdf">
        Dantzig, Fulkerson and Johnson</a> in 1954.
      </p>

      <p>
        In this example, we consider a salesman traveling in the US.
        The salesman starts in New York and has to visit a set of cities
        on a business trip before returning home. The problem then
        consists of finding the shortest tour which visits every city
        on the itinerary.
      </p>

      <p>
        We will formulate this problem as an integer program and implement
        it in Gurobi. The implementation will also demonstrate the use of
        <em>lazy constraints</em> in Gurobi.
      </p>
    </div>
    <div class="example_section" id="model">
      <h2><a name="model">Mathematical Model</a></h2>

      <p>
        We let the $n$ selected cities in the salesman's tour be the
        set of vertices $V$ of a graph. The set of edges $E$ of the
        graph corresponds to the different connections between each
        city. Since we can travel from any city to another, the graph
        is complete. That is, there is an edge between every pair of
        nodes. For each edge in the the graph we associate a binary
        variable
        \[
        x_{ij} = \left\{\begin{array}{ll}
               1 & \text{if edge $(i,j) \in E$ is in tour }\\
               0 & \mathrm{otherwise.}
              \end{array}\right.
        \]
        Since the edges are undirected, we have that $x_{ij} = x_{ji}$,
        and it suffices to only include edges with $i &lt; j$ in the model.
      </p>

      <p>
        We want to minimize the total distance travelled during the
        tour.  Therefore, we calculate the distance $d_{ij}$ between
        each pair of nodes $i$ and $j$. The total distance travelled
        is then the sum of the distances of the edges included in the
        tour
        \[
        \text{distance travelled} = \sum_{(i,j) \in E} d_{ij} x_{ij}.
        \]
      </p>

      <p>
        The tour should only pass through each city once. Therefore,
        each node in the graph should have exactly one incoming
        edge and one outgoing edge. In other words, for every node $i$
        exactly two of the $x_{ij}$ binary variables should be equal
        to 1. We write this constraint as

        \[
        \sum_{j \in V} x_{ij} = 2 \quad \forall i \in V.
        \]

        This constraint means that the saleman should enter
        and leave each city exactly once.
      </p>

      <p>
        With just this constraint on the number of edges in the tour
        entering and exit each node, we may produce solutions that are
        not connected tours. For example, the figure below shows a graph
        with 6 nodes and two disconnected subtours. The first subtour
        goes through nodes 1, 2, and 3, while the second subtour goes
        through nodes 4, 5, 6, and 7. Note that in this solution each node
        has exactly two edges in the tour incident to it. But there is
        no path for the salesman to travel between subtour. So we must
        add extra constraints to our model to eliminate these
        solutions.
      </p>

      <div id="subtourdemo"></div>

      <script src="https://cdnjs.cloudflare.com/ajax/libs/d3/3.5.5/d3.min.js"></script>
      <script>
        var width = 800,
            height = 250,
            size = width/6;

        var svgdemo = d3.select("#subtourdemo").append("svg")
                    .attr("width", width)
                    .attr("height", height);

        var lineG = svgdemo.append("g");

        var pointsG = svgdemo.append("g");

        var textG = svgdemo.append("g");

        var pointData = [[.1*width, .8*height], [.25*width, .2*height], [.35*width, .75*height],
                         [.5*width, .18*height], [.68*width, .22*height], [.72*width, .8*height],
                         [.45*width, .90*height]];

        var cycleData = [[0,1], [1,2], [2,0], [3,4], [4,5], [5,6], [6,3]];

        var points = pointsG.selectAll("circle")
                            .data(pointData)
                            .enter()
                            .append("circle")
                            .attr("cx", function(d) { return d[0]; })
                            .attr("cy", function(d) { return d[1]; })
                            .attr("r", 5)
                            .attr("stroke", "rgb(100, 0, 0)")
                            .attr("fill", "rgb(200, 0, 0)");

        var texts = textG.selectAll("text")
                           .data(pointData)
                           .enter()
                           .append("text")
                           .text(function(d,i) { return String(i+1); })
                           .attr("x", function(d) { return d[0] - 20; })
                           .attr("y", function(d) { return d[1] - 5; })
                           .attr("font-family", "sans-serif")
                           .attr("font-size", "18px")
                           .attr("fill", "black")
                           .attr("font-weight", 400)
                           .attr("text-anchor", "middle");

        var cycles = lineG.selectAll("line")
                          .data(cycleData)
                          .enter()
                          .append("line")
                          .attr("x1", function(d) { var j = d[0]; return pointData[j][0]; })
                          .attr("x2", function(d) { var j = d[1]; return pointData[j][0]; })
                          .attr("y1", function(d) { var j = d[0]; return pointData[j][1]; })
                          .attr("y2", function(d) { var j = d[1]; return pointData[j][1]; })
                          .attr("stroke", "black")
                          .attr("stroke-width", 2);
      </script>

      <p>
        To eliminate the subtours we add the following constraints:
        \[
        \sum_{i,j \in S, \, i \neq j} x_{ij} \leq \left\vert{S}\right\vert - 1, \quad \forall S \subset V, S \ne \emptyset
        \]
        These constraints require that for
        each proper (nonempty) subset $S$ of the set of cities $V$,
        the number of edges between the nodes of $S$ must be at most
        $\left\vert{S}\right\vert - 1$.
      </p>

      <p>
        Indeed if the number of edges were equal to $S$, it would be
        possible to form a subtour. For example, in the figure above
        the subset of nodes $S=\{1,2,3\}$ has 3 edges in the tour: $x_{13} = 1,
        x_{12} = 1, x_{23} = 1$. So
        \[
        \sum_{i,j \in \{1,2,3\}, i \ne j} x_{ij} =  3 > 2  = |\{1, 2, 3\}| - 1
        \]
        Thus, the subtour elimination constraint above is violated.
      </p>

      <p>
        So finally the integer program formulation becomes
        \[
        \begin{array}{ll}
        \text{minimize} & {\displaystyle \sum_{(i,j) \in E} d_{ij} x_{ij}} \\
        \text{subject to} & {\displaystyle \sum_{j \in V} x_{ij}} = 2 \quad \forall i \in V \\
                          & {\displaystyle \sum_{i,j \in S, \, i \neq j} x_{ij}} \leq  \left\vert{S}\right\vert - 1 \quad \forall S \subset V, S \ne \emptyset  \\
                          & x_{ij} \in \{ 0, 1 \}
        \end{array}
        \]
      </p>

      <p>If the set of cities $V$ is of size $n$, there are $2^n - 2$
        subsets $S$ of $V$, excluding $S = V$ and $S = \emptyset$. Instead
        of explicitally including a constraint
        \[
        \sum_{i,j \in S, \, i \neq j} x_{ij} \leq  \left\vert{S}\right\vert - 1
        \]
        for each $S$. We include above constraints implicitly as <em>lazy constraints</em>.
        That is, we generate and add these constraints to our model in a lazy fashion.
      </p>

      <p> Initially we will have no subtour elimination constraints in our model.
        When Gurobi finds a feasible solution that satisifies the other constraints,
        we compute the shortest cycle in edges included in the tour so far. If
        this cycle is of length $n$, then the model is solved. Otherwise, a
        cycle of length $m &lt; n$ defines a subtour, and a set $S$ with $|S| = m$.
        We then add the corresponding constraint to eliminate this subtour as
        a lazy constraint, and Gurobi continues solving this new modified model.
        This process continues until the shortest cycle is of length $n$, implying
        that all subtour elimination constraints have been satisfied.
      </p>

    </div>

    <div class="example_section" id="implementation">
      <h2><a name="implementation">Implementation</a></h2>
      <p>Below is the full implementation of the model (and the associated data) in
        Gurobi's Python interface:
      </p>
<examplecode>
import math
import random
from gurobipy import *


# Callback - use lazy constraints to eliminate sub-tours

def subtourelim(model, where):
  if where == GRB.callback.MIPSOL:
    selected = []
    # make a list of edges selected in the solution
    for i in range(n):
      sol = model.cbGetSolution([model._vars[i,j] for j in range(n)])
      selected += [(i,j) for j in range(n) if sol[j] &gt; 0.5]
    # find the shortest cycle in the selected edge list
    tour = subtour(selected)
    if len(tour) &lt; n:
      # add a subtour elimination constraint
      expr = 0
      for i in range(len(tour)):
        for j in range(i+1, len(tour)):
          expr += model._vars[tour[i], tour[j]]
      model.cbLazy(expr &lt;= len(tour)-1)


# Euclidean distance between two points

def distance(points, i, j):
  dx = points[i][0] - points[j][0]
  dy = points[i][1] - points[j][1]
  return math.sqrt(dx*dx + dy*dy)


# Given a list of edges, finds the shortest subtour

def subtour(edges):
  visited = [False]*n
  cycles = []
  lengths = []
  selected = [[] for i in range(n)]
  for x,y in edges:
    selected[x].append(y)
  while True:
    current = visited.index(False)
    thiscycle = [current]
    while True:
      visited[current] = True
      neighbors = [x for x in selected[current] if not visited[x]]
      if len(neighbors) == 0:
        break
      current = neighbors[0]
      thiscycle.append(current)
    cycles.append(thiscycle)
    lengths.append(len(thiscycle))
    if sum(lengths) == n:
      break
  return cycles[lengths.index(min(lengths))]

n = 50

# Create n random points

random.seed(1)
points = []
for i in range(n):
  points.append((random.randint(0,100),random.randint(0,100)))

m = Model()


# Create variables

vars = {}
for i in range(n):
   for j in range(i+1):
     vars[i,j] = m.addVar(obj=distance(points, i, j), vtype=GRB.BINARY,
                          name='e'+str(i)+'_'+str(j))
     vars[j,i] = vars[i,j]
   m.update()


# Add degree-2 constraint, and forbid loops

for i in range(n):
  m.addConstr(quicksum(vars[i,j] for j in range(n)) == 2)
  vars[i,i].ub = 0

m.update()


# Optimize model

m._vars = vars
m.params.LazyConstraints = 1
m.optimize(subtourelim)

solution = m.getAttr('x', vars)
selected = [(i,j) for i in range(n) for j in range(n) if solution[i,j] &gt; 0.5]
assert len(subtour(selected)) == n
</examplecode>
    </div>
    <div class="example_section" id="demo">
      <h2><a name="demo">Live Demo</a></h2>
      <p>
        The cities which are part of the tour have been highlighted: New York,
        Houston, San Francisco, Seattle, Minneapolis and Denver. Click "Compute
        Tour" to find the optimal tour using Gurobi.
      </p>
      <p>
        You can hover over cities to show their names and click to add and remove
        them from the tour.
      </p>
      <div id="demoarea">
      </div>
      <button class="pure-button" onclick="compute()">Compute Tour</button>
    </div>

    <p>
      <button class="pure-button" onclick="toggle_div()">Gurobi Log</button>
    </p>


    <script src="https://cdnjs.cloudflare.com/ajax/libs/d3/3.5.5/d3.min.js"></script>
    <script>
      function toggle_div() {
        var logfile = d3.select('#logfile');
        if (logfile.style("display") === "none") {
          logfile.style("display", "block");
        } else {
          logfile.style("display", "none");
        }

      }
    </script>

    <examplecode id=logfile>
    </examplecode>

    <div class="example_section" id="try">
      <h2><a href="#try" name="try">Try Gurobi for Free</a></h2>
      <p> We hope this example has taught you a bit about the
        traveling salesman problem and using Gurobi. We encourage you
        to try the example out for yourself with Gurobi.  To enable
        this, we provide easy access to a full-featured evaluation
        version of Gurobi.
      </p>
      <div class="col_5 column">
        <a href="http://www.gurobi.com/downloads/evaluation-request">
          <button class="red stack-button">
            <i class="fa fa-lg fa-line-chart"></i>
            Commercial Users: Free Evaluation Version
          </button>
        </a>
      </div>
      <div class="col_5 column">
        <a href="http://www.gurobi.com/downloads/download-center">
          <button class="red stack-button">
            <i class="fa fa-lg fa-line-chart"></i>
            Academic Users: Free Academic Version
          </button>
        </a>
      </div>
      <p>
        We are always happy to discuss your needs and answer your questions.
        Just <a href="http://www.gurobi.com/company/contact-us">contact us</a>
        at your convenience.
      </p>
    </div>

    <div style="min-height:100px"></div>
<!--[if gt IE 8]><!--><script src="//ajax.googleapis.com/ajax/libs/jquery/2.1.0/jquery.min.js"></script><!--<![endif]-->
<script src="jquery.nav.js"></script>
<script>
  $(document).ready(function() {
  console.log('calling onePageNav');
  $('#nav').onePageNav({scrollOffset:120});
  });
</script>
<script src="https://cdnjs.cloudflare.com/ajax/libs/d3/3.5.5/d3.min.js"></script>
<script src="http://d3js.org/topojson.v1.min.js"></script>
<script src="https://cdnjs.cloudflare.com/ajax/libs/spin.js/1.2.7/spin.min.js"></script>
<script>

// Hide Log File intially
d3.select('#logfile').style("display", "none");

var width = 960,
    height = 600;
    padding = 10;

var svg = d3.select("#demoarea").append("svg")
    .attr("width", width)
    .attr("height", height);

var voronoi = d3.geom.voronoi().clipExtent([[0, 0], [width, height]]);

// G objects for graphics
var mapG = svg.append("g");
var labelG = svg.append("g");
var lineG = svg.append("g");
var pointsG = svg.append("g");
var selectG = svg.append("g");
var voronoiG = svg.append("g");

// Data structures needed for airport data
var airports = [];
var airportInfo = [];
var airportIds = [];
var airportNum = 0;
var voroPath;

// Add the map
var projection = d3.geo.albers()
    .center([0, 39.7])
    .rotate([104.9, 0])
    .parallels([50, 60])
    .scale(1200)
    .translate([2*width / 5 , 0.5*height]);

var path = d3.geo.path()
    .projection(projection)
    .pointRadius(2);

// Create color gradient for the map
var mapColors = [];
var maxRed = 222, minRed = 49;
var maxGreen = 235, minGreen = 130;
var maxBlue = 247, minBlue = 189;
var numColors = 10;
for (var i = 0; i < numColors; i++) {
  var red = String( Math.round( minRed + (maxRed - minRed)*i/numColors ) );
  var green = String( Math.round( minGreen + (maxGreen - minGreen)*i/numColors) );
  var blue = String( Math.round( minBlue + (maxBlue - minBlue)*i/numColors ) );
  mapColors.push( "rgb(" + red + "," + green + "," + blue + ")");
}

d3.json("US_TSP.json", function(error, us) {
  var states = topojson.feature(us, us.objects.states);
  var places = topojson.feature(us, us.objects.placesUSA5);

  mapG.selectAll(".subunit")
      .data(states.features)
      .enter().append("path")
      .attr("fill", function(d,i) { return mapColors[i % 9];})
      .attr("stroke", "white")
      .attr("stroke-width", 1)
      .attr("opacity", .8)
      .attr("d", path);

  mapG.append("path")
      .datum(places)
      .attr("d", path)
      .attr("class", "place");

  for (var i = 0; i < places.features.length; i++) {
    airports.push(projection(places.features[i].geometry.coordinates));
    airportInfo.push(places.features[i].properties.name);
  }
  console.log('number of places', places.features.length)
  console.log('airports', airports);

  // Filter because some of the airports are also in Canada, Mexico etc...
  var dataAirports = voronoi(airports).filter(nonzeroArray);
  console.log('data', dataAirports);
  console.log('airIds', airportIds);

  labelG.selectAll(".place-label")
      .data(places.features)
      .enter().append("text")
      .attr("class", "place-label")
      .attr("transform", function(d,i) { return "translate(" + airports[i] + ")"; })
      .attr("x", function(d,i) { return airports[i][0] > -1 ? 6 : -6; })
      .attr("dy", ".35em")
      .style("text-anchor", function(d,i) { return airports[i][0] > -1 ? "start" : "end"; })
      .text(function(d) { return d.properties.name; })
      .attr("id", function(d,i) { if (airportIds[i] > -1) { return "p" + airportIds[i]; } else {return "n" + i;} })
      .attr("opacity", 0);

  var voroPath = voronoiG.selectAll("polygon")
                        .data(dataAirports)
                        .enter()
                        .append("polygon")
                        .attr("points", function(d) { return d; })
                        .attr("fill", "black")
                        .attr("opacity", 0)
                        .attr("id", function(d,i) { return i;})
                        .attr("stroke", "black")
                        .attr("stroke-width", 2)
                        .attr("stroke-opacity", 1)
                        .on("mouseover", mouseover)
                        .on("mouseout", mouseout)
                        .on("mousedown", mousedown);
});

function nonzeroArray(value) {
  if (value.length > 0) {
    airportIds.push(airportNum);
    airportNum++;
  } else {
    airportIds.push(-1);
  }
  return value.length > 0;
}

function indexArray(array, point) {
  var p0 = point[0];
  var p1 = point[1];
  var index = -1;
  for (i = 0; i < array.length; i++) {
    pointprime = array[i];
    if (p0 === pointprime[0] && p1 === pointprime[1]) {
      index = i;
      break;
    }
  }
  return index;
}

function polygon(d) {
  return "M" + d.join("L") + "Z";
}

function mouseover() {
  var id = d3.select(this).attr("id");
  var place = "#p" + id;
  d3.select(place).attr("opacity", 1);
}

function mouseout() {
  var id = d3.select(this).attr("id");
  var place = "#p" + id;
  d3.select(place).attr("opacity", 0);
}

function mousedown() {
  lineG.selectAll("line").remove("line");

  var id = d3.select(this).attr("id");
  id = parseInt(id);
  console.log('id', id);
  var num = airportIds.indexOf(id);
  console.log('num',num);
  var point = airports[num];
  console.log(point);

  if (indexArray(nodes, point) > -1) {
    console.log('point already exists!');
    nodes.splice(indexArray(nodes, point),1);
  } else {
    nodes.push(point);
  }

  // Draw the nodes
  pointsG.selectAll("circle").remove("circle");

  var circles = pointsG.selectAll("circle")
                 .data(nodes)
                 .enter()
                 .append("circle");

  circles.attr("cx", function(d, i) {
            return nodes[i][0];
            })
            .attr("cy", function(d, i) {
            return nodes[i][1];
            })
            .attr("r", 4)
            .attr("fill",  "steel blue")
            .attr("stroke", "steel blue")
            .attr("stroke-opacity", 0.5)
            .attr("stroke-width", 3);
}

// New York, Houston, SF, Seattle, Minneapolis, Denver coordinates
var nodes = [[868.5341864421549, 170.05210236092614], [568.2933565139988, 485.90016899260206],
             [88.19112235025932, 302.28594437236313], [139.08656500094133, 106.64803640280218],
             [557.3170513891673, 176.76299706303132], [382.60278181309786, 299.11173993840816]];

// Add initial cities in tour
var circles = pointsG.selectAll("circle")
                 .data(nodes)
                 .enter()
                 .append("circle");

  circles.attr("cx", function(d, i) {
            return nodes[i][0];
            })
            .attr("cy", function(d, i) {
            return nodes[i][1];
            })
            .attr("r", 4)
            .attr("fill",  "steel blue")
            .attr("stroke", "steel blue")
            .attr("stroke-opacity", 0.5)
            .attr("stroke-width", 3);

function compute() {
  console.log('vertices', nodes);
  if (nodes.length < 3) {
    alert("At least 3 cities are required to compute a tour. Add more!")
    return;
  }
  d3.json('airlineTSP.py')
    .header('Content-Type', 'application/json')
    .post(JSON.stringify({'vertices': nodes}), serverResponse);
}


function serverResponse(error, data) {
   console.log('serverResponse');
   console.log('data', data);
   if (!error) {
      if ('tour' in data) {
        var solution = data.tour;
        if (solution[0] === "error") {
          alert("10 second solve time limit exceeded.");
          return;
        }
        var tour = solution[0];
        console.log('tour', tour);

        var logMsg = solution[1]; // Log message to display

        d3.select('#logfile').html(logMsg);

        // Get list of edges of tour
        var src = tour[0];
        links = [];
        for (var i = 1; i < tour.length; i++) {
          links.push([src, tour[i]]);
          src = tour[i];
        }
        links.push([tour[tour.length-1], tour[0]]);

        console.log('links', links);

        // Draw lines between points
        lineG.selectAll("line").remove("line");

        var lines = lineG.selectAll("line")
                         .data(links)
                         .enter()
                         .append("line");

        lines.style("stroke", "purple")
             .style("opacity", 0.5)
             .style("stroke-width", 3)
              .attr("x1", function(d,i) {
                var node1 = nodes[links[i][0]];
                var x1 = node1[0]
                return x1;
              })
              .attr("y1", function(d,i) {
                var node1 = nodes[links[i][0]];
                var y1 = node1[1]
                return y1;
              })
              .attr("x2", function(d,i) {
                var node2 = nodes[links[i][1]];
                var x2 = node2[0]
                return x2;
              })
              .attr("y2", function(d,i) {
                var node2 = nodes[links[i][1]];
                var y2 = node2[1]
                return y2;
              })
      }
   }
}

</script>
